Characterizing the dispersion behavior of poly-atomic magnetic metamaterials

The propagation of magnetoinductive (MI) waves across magnetic metamaterials known as magnetoinductive waveguides (MIWs) has been an area of interest for many applications due to the flexible design and low-loss performance in challenging radio-frequency (RF) environments. Thus far, the dispersion behavior of MIWs has been limited to mono- and diatomic geometries. In this work, we present a recursive method to generate the dispersion equation for a general poly-atomic MIW. This recursive method greatly simplifies the ability to create closed-form dispersion equations for unique poly-atomic MIW geometries versus the previous method. To demonstrate, four MIW geometries that have been selected for their unique symmetries are analyzed using the recursive method. Using applicable simplifications brought on by the geometric symmetries, a closed-form dispersion equation is reported for each case. The equations are then validated numerically and show excellent agreement in all four cases. This work simultaneously aids in the further development of MIW theory and offers new avenues for MIW design in the presented dispersion equations.

As previously discussed, the theoretical analysis of one-dimensional MIWs has thus far been limited to monoand di-atomic MIWs 9,13 .The expansion from the one-to two-element case allowed for expanded bandwidths, two separable bands of operation, reduced path loss, and the supporting of both forward and backward waves.These benefits came at the cost of increased difficulty in impedance matching due to the number of elements that must be accounted for when synthesizing the matching impedance 13 and an increased analysis complexity.In particular, the difficulty in transducer design for the di-atomic case has been the primary reason that they remain mostly unused in MIW applications.Despite the even further increased challenges in impedance matching, expanding to a general poly-atomic case (K resonant elements per periodic unit) is still of interest.The expansion would allow for up to K separable bands of operation and potential for even further reduced path loss and expanded bandwidths which would benefit many of the applications of MIWs described previously.One reason this expansion has yet to be realized in theory, besides the foreseen challenges in implementation and impedance matching, are the computational challenges facing the generation of the dispersion equations.The traditional analysis method relies on the determinant of the impedance matrix which grows in size as the square of the number of elements in the periodic unit (e.g., K = 3 has 9 entries, K = 4 has 16 entries, etc.) which becomes computationally intensive for large values of K.An alternative approach to this method for specific geometries is to examine the geometry as an infinite two-dimensional array as described by Shamonina 9 and assume a standing wave solution in one direction.This method requires a strong working knowledge of 2D MIW arrays and also breaks down in the general case.So, while there exists interest in expanding the theory beyond the di-atomic case due to the potential for improved performance, it has not yet been explored due to the cumbersome nature of the traditional analytical method and the inability to generalize the truncated array method.
In this work, we create a novel recursive method to generate a dispersion equation for a general poly-atomic (K resonant elements per periodic unit) MIW through the use of the nearest neighbor approximation and coupled chain analysis.This method greatly simplifies the generation of new closed-form dispersion equations compared to traditional methods.Using this new recursive method, we generate closed-form dispersion equations for four geometries with unique symmetries and applicable simplifications, and then validate these equations numerically as both a proof-of-concept for the utility of the new method and an expansion of the current MIW design space.

Derivation of recursive method
Figure 1a shows a truncated diagram of a general poly-atomic MIW.In this general scenario, we assume that each resonant element is identical (as is often the case), the MIW is periodic in the direction of propagation (i.e., the n − 1 th unit cell is identical to the nth unit cell and so forth), and that the MIW is infinite in length in the direction of propagation such that terminal reflections can be ignored.Additionally, we assume the coupling between non-nearest neighbor periodic units is significantly smaller than nearest neighbor coupling and can, thus, be ignored.Note that this is approximation is used only to simplify the presentation of the analytical work.We will show later on that expanding the method to include non-nearest neighbor interactions between periodic units is simple and follows from previous theory.As seen, the elements that make up each unit can be rotated at any angle (in 3D space) and can be shifted by up to +/− half of the geometric period of the MIW relative to the first element which we select as the 0 reference point.Table 1 contains the dictionary for the various terms used throughout the derivation.
The proposed derivation of the relationship between the complex propagation constant, γ , and the radial frequency, ω , for a poly-atomic MIW follows a generalized form of the coupled chain analysis conducted for the two-element case 13 .We begin by invoking Kirchoff 's Voltage Law on each element in the nth unit, assuming that the non-nearest neighbor coupling in the direction of propagation is negligible.This generates three types of expressions due to the finite number of elements: (1) an expression for the 1st element, (2) an expression for any kth element, and (3) an expression for the Kth element.For illustrative purposes, we will focus just on the 1st element, with all others following a similar procedure and form.Assuming no external excitation on the nth loop: We then look for wave solutions of the form I n,k = I k e −njγ p which corresponds to a traveling wave in the positive n direction with γ = β − jα .After simplification, we reach the following form for the first element: Similar expressions are created for the kth and Kth elements as well.
From these expressions, we convert the system of equations to matrix form such that V K = Z K I K .Here, V K is a K × 1 0 vector and I K contains the K current constants noted previously.To simplify the notation of the impedance matrix, we set a i = Z + 2Z i cos(γ p) and b ik = Z (n,k),(n,i) + e jγ p Z (n−1,k),(n,i) + e −jγ p Z (n+1,k),(n,i) .Thus the impedance matrix for a general poly-atomic MIW is: (1) (n,1) + e jγ p Z (n−1,i),(n,1) + e −jγ p Z (n+1,i),(n,1) www.nature.com/scientificreports/To generate non-trivial solutions to the equation V K = Z K I K , we require that |Z K | = 0 .This condition relates the propagation constant ( γ ) and the radial frequency ( ω ), thus creating a dispersion equation for any MIW.The problem with generalizing this condition is that each element in Z K is unique, i.e., generalizing the condition means generalizing the determinant.This is a problem that is already well studied and has several solutions.For example, we can use the Laplace expansion to generate the determinants through the cofactor matrix.As (3) Figure 1.Diagram of poly-atomic MIWs with K elements per unit cell.Each rectangle represents a resonant element in a specific orientation.The elements are labeled by their element number in the periodic unit (k) and their unit number (n).Each element is identical and the direction of propagation is horizontal.(a) General case with all elements that share an element number having the same orientation as one another.(b) Case (1) with each element in the periodic unit aligned in orientation and position and with uniform spacing between elements within the periodic unit.(c) Case 2) where each element in the periodic unit is aligned in orientation and the spacing between elements within the periodic unit is uniform.An alternating shift between elements in the periodic unit is introduced.(d) Case (3) has each element in the periodic unit aligned in position and has the uniform spacing between elements in the periodic unit, but an alternating 90 degree rotation of elements is introduced.(e) Case (4) has uniform spacing within the periodic unit but an alternating rotation and shift of elements.
previously discussed, this is the typical approach used to generate closed-form dispersion equations for any MIW 9,13 .Our goal is to create a simplified method of generating these dispersion equations and to do so, we must enforce further simplifications.
To simplify the expression, we further apply the nearest neighbor approximation such that we only account for nearest neighbor coupling in the direction of propagation and within each periodic unit.Physically, this means that the elements within our periodic unit are weakly coupled to one another such that the coupling between elements (n, k) and (n, k ± 1) is much larger than the coupling between elements (n, k) and (n, k ± 2) .Instead of working through the derivation process again, we can simply see that this approximation leads to all b ik = 0 for |i − k| > 1 .This leads to a new impedance matrix: With many 0 terms, there is now potential to further generalize the determinant condition beyond the typical approaches.Through the use of the Laplace expansion, the first few dispersion equations are: From Eqs. ( 5)-( 7), a clear pattern emerges that does continue to higher values of K.For a general poly-atomic MIW, we can now determine the dispersion equation recursively as: Where |Z K | is the determinant of matrix Z K .Although this recursive method will generate a general disper- sion equation for any value of K, we can also use this method to examine simplified geometries to generate a closed-form expression for their dispersion equations.We will examine four (4) distinct cases: Case (1) Fully Uniform Array, Case (2) Uniform Array with Alternating Shift, Case (3) Uniform Array with Alternating Rotation, Case (4) Uniform Array with Alternating Rotation and Shift.While we will maintain the nearest neighbor approximation between periodic units for each Case to simplify the analysis, this recursive method can still be used in the case of significant non-nearest neighbor coupling between periodic units with an extension of the a i and b ik definitions.( 4) Dictionary of mathematical symbols and their corresponding meaning for the derivation and application of the recursive method to generate closed-form dispersion equations.

Symbol
Parameter Mutual inductance between elements (n,k) and (u,q) Mutual impedance between elements (n,k) and (u,q) I n,k Current on loop (n,k) www.nature.com/scientificreports/

Application
For Case (1), we take the restriction to the extreme and enforce rotation and position alignment and uniform spacing of all elements as shown in Fig. 1b.Mathematically, this leads to . In terms of the matrix terms, we now have ),(n,1) + 2Z (n±1,2),(n,1) cos(γ p) .With these new expressions, we again look at the first few dispersion equations: With Eqs. ( 9)-( 14) completed, a pattern emerges as a result of the symmetry.Based on this, the closed-form dispersion equation for a poly-atomic MIW with identical, aligned elements and uniform spacing within the periodic unit is: From this expression, the polynomial nature of the dispersion equations for MIWs becomes very clear.The introduction of additional elements in each periodic unit leads to additional solutions to the polynomial expression.
In turn, this leads to additional propagation constants such that a poly-atomic MIW may have up to K separable frequency bands of operation.For Case (2), we generalize the expression further by relaxing some of the uniformity requirements.First, we allow for a periodic shift of every other element in the unit as shown in Fig. 1c.We still have but we have different symmetries for the diagonal terms depending on whether k is even or odd.For even k, we have: For odd k, the symmetry dictates: T h i s g i v e s u s t w o o f f -d i a g o n a l t e r m s a n d o n e o n -d i a g o n a l t e r m : ,(n,1) + e jγ p Z (n+1,2),(n,1) + e −jγ p Z (n−1,2),(n,1) .Following the same process with the dispersion equa- tions, a similar pattern emerges which gives the following general dispersion equation: The equation is very similar to Eq. ( 15) as anticipated and also simplifies to Eq. ( 15) when the geometric shift is 0 or ±p/2 , such that b 1 = b 2 .
Case (3) is the dual of Case (2) in that we maintain the off-diagonal symmetry while introducing periodicity in the on-diagonal entries from the initial case.This can be easily accomplished through an alternating 90-degree rotation between aligned elements in each unit, as shown in Fig. 1d.By rotating the elements, we now have an additional Z 1 = Z 3 = Z 5 . . .and Z 2 = Z 4 = Z 6 = . . . .This implies the two on-diagonal terms a 1 = Z + 2Z 1 cos(γ p) and a 2 = Z + 2Z 2 cos(γ p) .Examining Eq. ( 8), we see that we now have sepa- rate expressions when K is even or odd.Because of the center alignment of the elements, we still maintain b ik = b ki = b = Z (n,2),(n,1) + 2Z (n±1,2),(n,1) cos(γ p) .The new recursive relation is: www.nature.com/scientificreports/While the even/odd nature of the relationship complicates the analysis, the same procedure can be followed to generate a closed-form solution that is similar to both Eqs. ( 15) and ( 20): It is important to note that by introducing a rotated element into the unit cell, there is now the possibility of both forward and backward waves propagating on the structure.This phenomenon has been explored in depth in previous literature 13 .Finally, in Case (4), we combine both Case ( 2) and ( 3) by introducing a periodic shift and 90-degree rotation within the unit cell, as shown in Fig. 1e.Without further analysis, it is clear that the symmetries will include a combination of the previous two cases such that are the only terms in the impedance matrix.Following the derivations of Eqs. ( 20) and (23), this leads to a dispersion equation of the form: As previously mentioned, the presented dispersion equations in Eqs. ( 15), ( 20), (23), and (28) rely on the nearest neighbor approximation in all directions for each resonant element.While we cannot relax the approximation within the unit cell while maintaining the form of the expressions, we can relax the approximation between unit cells with relative ease.As an example, when considering U nearest neighbors between unit cells of the uniform and aligned MIW shown in Fig. 1b, the dispersion equation remains the same, just with new definitions of a and b.This foundation can be easily applied for each of the other cases explored throughout this work.

Validation
To validate each analytical equation, we use the CST Studio finite-element method, frequency domain, 3D fullwave electromagnetic solver 49 .This method was selected to demonstrate that the relatively simple dispersion equations can predict MIW performance with a high-level of accuracy despite not accounting for effects such as non-nearest neighbor coupling, electric coupling, and radiation that are present in both the full-wave simulation and in reality.While experimental validation is outside the scope of this work, we note strong agreement between similar simulations and experiments have been previously demonstrated 27,28,30 .As such, the presented simulation results are expected to be a good approximation of experimental results.
For each MIW in each Case, the resonant elements are 9.1 × 3.5 cm copper rectangular loops made of 30 American Wire Gauge wire and loaded with 56 pF capacitors to achieve resonance at 41.3 MHz.This design corresponds to equivalent circuit parameters of R = 0.92 , L = 260.3nH, and C = 56 pF.The physical gap between units in all scenarios is set at 0.25 cm and the loops are oriented in such a way that the overall geometric period is 3.75 cm.The distance between elements is set to 0.75 cm.Both of these distances were selected to reduce the effect of non-nearest neighbor coupling to meet our base assumptions for each equation.Table 2 contains the appropriate mutual inductance values for each Case.
To simulate each Case in CST Studio, we truncate the MIWs to 11 units, place the structure in free space, and excite the first element in the first unit with a 0.5 W sinusoidal signal.S-Parameters are then generated using a reference impedance of 50 .Because of the selected simulation method and lack of a mode-selective and wellmatched transducer for poly-atomic MIWs, the wave modes cannot be easily separated.As such, the S-Parameters are presented separately from the calculated attenuation and phase constants.Instead of directly comparing the analytical model results to equivalent propagation constants calculated from the S-Parameters, the bandwidth and loss predicted by the analytical model will be compared to the performance shown in the corresponding (n,1) + e jγ p Z (n−1,2),(n,1) + e −jγ p Z (n+1,2),(n,1) ,(n,1) + e jγ p Z (n+1,2),(n,1) + e −jγ p Z (n−1,2),(n,1) Case (1), as shown in Fig. 1b with the dispersion equation defined in Eq. ( 15), is analyzed for the K = 3 and K = 4 cases in Fig. 2a,b.Looking at the attenuation constants ( α ), we first see the expected K separate solutions from the analytical model.For both K = 3 and K = 4 , the first solution corresponds to the lowest attenuation across the examined frequency range.Figure 3a For the phase constant ( β ), the bandwidth of each solution is defined by the frequencies where the deriva- tive of the phase constant with respect to frequency is nonzero (i.e., exhibits a non-zero group velocity).For the S-Parameter results, the bandwidth definition is not as clear.We will define the number of bands and their corresponding bandwidths with the following procedure: (1) determine the value and frequency of the minimum attenuation, (2) from this point, increase and decrease the frequency until the attenuation is 20 dB less than the minimum attenuation found in (1) and set these two frequency points as the high and low-frequency cutoffs for the current band.(3) repeat steps (1) and ( 2), removing the previous band from the analysis.This process repeats a maximum of K times, corresponding to the K possible solutions.In the K = 3 scenario for Case (1), the analytical results predict three passbands with the corresponding fractional bandwidths (with respect to the resonant frequency): 0.27, 0.30, and 0.30.Following our procedure, the S-Parameter results indicate fractional bandwidths of 0.22, 0.22, and 0.10 for each respective passband.While these bandwidths are significantly smaller than the analytical results, the analytical values do not account for the significant overlap in frequency between each solution.Instead we determine the overall bandwidth of the structure regardless of mode for the analytical results and compare this to the overall bandwidth of the simulated results.The analytical model predicts an overall bandwidth of 0.64 while the simulated results have an overall bandwidth of 0.56.This slight decrease in bandwidth is expected as the -20 dB definition of bandwidth used for the simulated results is stricter than the definition used in the analytical results which allows for extremely high levels of loss at the edges of the bands.Increasing K to 4 leads to similar results with the analytical model predicting bandwidths of 0.30, 0.30, 0.30 and 0.27 for each band and an overall bandwidth of 0.71 while the simulated results have bandwidths of 0.18, 0.16, 0.08, and 0.17 with an overall bandwidth of 0.58.
For Case (2), we introduce a quarter period (0.9325 cm) shift in every other element in the periodic unit.The analytical results of this second case, as shown in Fig. 1c that has a dispersion equation defined by Eq. ( 20), are shown in Fig. 2c,d and the simulated results are in Fig. 3c,d.In the K = 3 scenario, only the first analytical passband is well represented by the simulated results while the upper bands are assimilated into a single passband due to large fluctuations in loss after the resonant frequency.This problem does not persist when K is increased to 4. Further examining the attenuation results, we see that the simulation upper frequency cutoff is decreased relative to the analytical model for both K = 3 and K = 4 .This difference is likely caused by stronger relative coupling between non-nearest neighbors within the periodic unit when compared to the first case.By shifting every other element, the coupling between the kth and k + 1 th elements is decreased while the coupling between the kth and k + 2 th elements is unchanged.This decreases the validity of our base assumption that non-nearest neighbor interactions are negligible.
For the phase constants, the K = 3 analytical model predicts passband with bandwidths 0.24, 0.30, 0.37 and an overall bandwidth of 0.58.Similarly, the K = 4 model predicts bandwidths of 0.27, 0.27, 0.34, 0.37 and an overall bandwidth of 0.64.The separation between bands in the simulated results is poor as anticipated by the significant overlap between solutions in the analytical model as such the individual passband bandwidths are not meaningful.The overall bandwidth for both K = 3 and K = 4 align well with the analytical results, at 0.49 and 0.51 respectively with the primary decrease in bandwidth relative to the model coming from the decrease in upper frequency cutoff as previously described.
For Case (3), we rotate every other element in the periodic unit by 90 degrees while maintaining uniformity in center alignment.This implies both positive and negative mutual coupling between the periodic units, which introduces both forward and backward waves as previously discussed.Figure 2e,f shows the results for this case based on the dispersion equation presented in Eq. ( 23) and Fig. 3e,f contain the simulated results.
We see strong agreement for the K = 3 and K = 4 attenuation constants, with a slight discrepancy in the center of the passband where both simulation results present a slight increase in attenuation rather than the expected smooth behavior shown in the analytical model results.This is likely caused by non-nearest neighbor coupling within the periodic unit being significantly stronger than the nearest neighbor coupling within the www.nature.com/scientificreports/periodic unit due to the 90 degree rotation.By placing the elements perpendicular to one another while maintaining their center alignment, we theoretically set the mutual coupling value between them to 0. As such, any non-zero coupling between non-nearest neighbor will have more impact on the results than other scenarios.
For the phase constant, we see the expected forward and backward waves in the analytical solutions.This is demonstrated by both negative and positive sloped curves across the phase constant.Regarding each passband, we have bandwidths of 0.24, 0.34, and 0.24 for the K = 3 case and bandwidths of 0.20, 0.27, 0.24, and 0.24 for the K = 4 case.As the phase constants lie essentially across identical frequencies, the overall predicted bandwidths of the structures are very close to the bandwidths of the individual solutions, at 0.34 for both cases.Likewise, the overall bandwidths in the simulated results are the same for both cases, at 0.28 which shows overall strong agreement in bandwidth between simulated and analytical results.
To model Case (4), we start at the geometry of Case (3) and shift the alternating rotated elements by a quarter period relative to the first element in the direction of propagation.This scenario is shown in Fig. 1e and modeled by Eq. (28).The analytical results are presented in Fig. 2g,h and simulated results in Fig. 3g,h.
As seen, the dispersive behavior is similar to Case (3) for both K = 3 and K = 4 .This is because of the relatively weak coupling between adjacent elements in the periodic unit that is maintained by only shifting the elements by a small distance.Comparing to the simulation, we also see a very similar level of agreement for  3) with alternating 90 degree rotations between aligned elements in the periodic unit with uniform spacing.K = 3 and K = 4 for (e) and (f).(g, h) Case (4) with alternating shifts and 90 degree rotations between elements in the periodic unit with uniform spacing.K = 3 for (g) and K = 4 for (h).both K = 3 and K = 4 .Again, the simulated attenuation constant is nearly identical to the analytical solutions, besides an increase in attenuation in the center of the passband.We can again explain this by the effects of nonnearest neighbor coupling.
For the phase constants, we again see strong agreement in the overall bandwidth between the simulated and analytical results for both K = 3 and K = 4 The model predicts an overall bandwidth of 0.37 for both K = 3 and K = 4 while the simulated results show overall bandwidths of 0.24 and 0.29 respectively.The larger difference in bandwidth for the K = 3 case is primarily caused by a significant decrease in upper frequency cutoff compared with the K = 4 case.

Conclusion
In this work, we have significantly expanded the MI wave guiding theory by defining a new recursive method of determining the dispersion equation for a general poly-atomic (K resonant elements per periodic unit) MIW.This recursive method relies on the validity of the nearest neighbor approximation within the periodic unit but can easily be expanded to include higher order coupling effects in the direction of propagation.Using this method, four cases were examined, each with unique symmetries and applicable simplifications.For each case, a closed-form dispersion equation was generated using the recursive method to greatly reduce the complexity of analysis.Each case was then validated through the use of simulation results for K = 3 and K = 4 .Results showed strong agreement between the analytical and numerical attenuation and phase behaviors over frequency.Discrepancies in results are likely caused by non-nearest neighbor coupling within the periodic unit that are not accounted for in the new recursive method and subsequently generated closed-form equations.
Through the new dispersion equations and recursive dispersion equation generation method, the MIW design space has been expanded significantly.In particular, the theory was previously limited to mono-or di-atomic MIWs but has now been expanded to any number of elements per period.The new closed-form equations expand the specific theory of MI wave propagation while the new method expands the tools available to those interested in specific and unique geometries for their applications such as imaging, communications, or sensing.While these new ideas rely on the validity of the nearest neighbor approximation within the periodic unit, they still serve as a useful approximation for cases where non-nearest neighbor coupling is significant, particularly for large value of K where the typical analysis becomes cumbersome.Future work will include validating the results presented here experimentally, expanding the analysis to include the synthesis of the required terminal impedance to reduce reflections in poly-atomic MIWs, and the design of mode-selective, well-matched transducers for di-atomic and poly-atomic MIWs.
,b show the corresponding S-Parameter results.The |S 21 | results clearly show the same K separated passbands with relatively large nulls separating the bands.Additionally, the loss increases per band as the frequency increases for each case in both the S-Parameter and attenuation constant results.

Figure 2 .
Figure 2.Frequency vs complex propagation constant results.Left image is imaginary attenuation constant ( α ) and the right is the real phase constant ( β ). (a, b) Case (1) with uniform spacing, orientation and position in the periodic unit with K = 3 and K = 4 respectively.(c, d) Case (2) K = 3 and K = 4 with uniform orientation and spacing but alternating shift in the periodic unit.(e, f) Case (3) with alternating 90 degree rotations between aligned elements in the periodic unit with uniform spacing.K = 3 and K = 4 for (e) and (f).(g, h) Case (4) with alternating shifts and 90 degree rotations between elements in the periodic unit with uniform spacing.K = 3 for (g) and K = 4 for (h).

Figure 3 .
Figure 3. S-Parameter versus normalized frequency results.|S 11 | results are shown with dashed lines and |S 21 | results are shown with solid lines.(a, b) Case (1) with uniform spacing, orientation and position in the periodic unit with K = 3 for for (a) and K = 4 for (b) respectively.(c, d) Case (2) K = 3 and K = 4 with uniform orientation and spacing but alternating shift in the periodic unit.(e, f) Case (3) with alternating 90 degree rotations between aligned elements in the periodic unit with uniform spacing.K = 3 and K = 4 for (e) and (f).(g, h) Case (4) with alternating shifts and 90 degree rotations between elements in the periodic unit with uniform spacing.K = 3 for (g) and K = 4 for (h).

Table 2 .
All mutual inductance values used in the validation of Case (1)-Case (4).